Data

The data we have are counts of deaths per episode of Game of Thrones.

GoT <- read_xlsx("GoT_Deaths.xlsx", col_names = T)
head(GoT)

## # A tibble: 6 x 3
##   Season Episode Count
##    <dbl>   <dbl> <dbl>
## 1      1       1     4
## 2      1       2     3
## 3      1       3     0
## 4      1       4     1
## 5      1       5     5
## 6      1       6     4

A quick histogram of the data:

ggplot(GoT, aes(x = Count))+
  geom_histogram(bins = 25, fill = "#deebf7") +
  labs(x = "Number of deaths", y = "Count")

Because we have count data, it is only natural to model the data with a Poisson distribution. If we take the empirical mean of the data as our estimate for \(\hat{\lambda}\), what does the sampling distribution look like?

y <- GoT$Count
n <- length(y)

# Finish: obtain empirical mean
ybar <- mean(y)
sim_dat <- rpois(n, ybar)
qplot(sim_dat, bins = 20, xlab = "Simulated number of deaths", fill = I("#9ecae1"))

# create data frame for side-by-side histograms
df <- rbind(data.frame(y, "observed") %>% rename(count = 1, data = 2),
            data.frame(sim_dat, "simulated") %>% rename(count = 1, data = 2))
ggplot(df, aes(x = count, fill = data))+
  geom_histogram(position = "dodge", bins = 20) +
  scale_fill_brewer()  +
  labs(x = "Number of deaths", y = "Count")

Plotting the histogram of the data next to the histogram of the sampled data using \(\hat{\lambda} = \bar{y}\), we see that there are many more 0-valued observations in the observed data than there are in the simulated data. We suspect that the Poisson model will not be a good fit for the data, but we will proceed anyway.

Fit Poisson model

We have learned that the Gamma distribution is conjugate for Poisson data. Here, we have chosen the prior \(\lambda \sim Gamma(10,2)\) based on prior knowledge that Game of Thrones is a show filled with death. We run the model in Stan. We save the posterior sampled values of \(\lambda\) in the variable `lambda_draws’.

stan_dat <- list(y = y, N = n)
fit <- stan("poisson-simple.stan", data = stan_dat, refresh = 0, chains = 2)

## Warning in readLines(file, warn = TRUE): incomplete final line found on '/
## Users/jb/Work/Teaching/STA601/labs/undergraduate/poisson-simple.stan'

lambda_draws <- as.matrix(fit, pars = "lambda")

We use the function mcmc_areas() to plot the smoothed posterior distribution for \(\lambda\). The line represents the posterior mean, and the shaded area represents 90% of the posterior density. You can play around with changing the `prob’ parameter in mcmc_areas(). We also report the posterior mean and 90% HPD interval using this Poisson-Gamma model.

mcmc_areas(lambda_draws, prob = 0.9)

print(fit, pars = "lambda")

## Inference for Stan model: poisson-simple.
## 2 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=2000.
## 
##        mean se_mean   sd 2.5%  25% 50%  75% 97.5% n_eff Rhat
## lambda    4    0.01 0.23 3.59 3.84   4 4.15  4.48   820    1
## 
## Samples were drawn using NUTS(diag_e) at Fri Jul 12 12:29:30 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

How does the posterior mean for \(\lambda\) compare to the sample mean?

Graphical Posterior Predictive Checks (PPC)

Let us generate posterior predictive samples using the posterior values of \(\lambda\).

y_rep <- matrix(NA, nrow = length(lambda_draws), ncol = n)
for(i in 1:length(lambda_draws)){
  # Finish generating posterior predictive samples
  y_rep[i,] <- rpois(n, lambda_draws[i])
}

We now step through some graphical posterior predictive checks. Here we plot the histogram of the observed counts to several of the posterior predictions. What do you notice?

ppc_hist(y, y_rep[1:8, ], binwidth = 1)

We can compare the density estimate of \(y\) to the density estimates for several (60) of the \(y_rep\)s. What do you notice here?

ppc_dens_overlay(y, y_rep[1:60, ])

In our first histogram of the data, we noticed the high number of 0 counts. Let us calculate the proportion of zeros in \(y\), and compare that to the propotion of zeros in all of the posterior predictive samples.

prop_zero <- function(x){
  mean(x == 0)
} 
prop_zero(y)

## [1] 0.1643836

ppc_stat(y, y_rep, stat = "prop_zero")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

We can plot the means and variances for all of the simulated datasets, along with the observed statistics.

ppc_stat_2d(y, y_rep, stat = c("mean", "var"))

ppc_error_hist(y, y_rep[1:4, ], binwidth = 1) + xlim(-15, 15)

## Warning: Removed 8 rows containing missing values (geom_bar).

Based on on these PPCs, does this model appear to be a good fit for the data?

Poisson Hurdle Model

How should we account for the high frequency of zeros? As it turns out, there is much literature on how to work with this sort of data. We can fit a zero-inflated Poisson model, where the data is modelled as a mixture of a Poisson distribution with a point mass at zero. Here, we fit a hurdle model: \(y_i = 0\) with probability \(\theta\), and \(y_i > 0\) with probability \(1-\theta\). Conditional on having observed \(y_i\) nonzero, we model \(y_i\) with a Poisson distribution which is truncated from below with a lower bound of 1. For the truncated Poisson, we use the same \(Gamma(10,2)\) prior as in the simple Poisson model above.

fit2 <- stan("poisson-hurdle.stan", data = stan_dat, chains = 2, refresh = 0)

## Warning in readLines(file, warn = TRUE): incomplete final line found on '/
## Users/jb/Work/Teaching/STA601/labs/undergraduate/poisson-hurdle.stan'

## Warning: There were 8 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded

## Warning: Examine the pairs() plot to diagnose sampling problems

print(fit2, pars = c("lambda", "theta"))

## Inference for Stan model: poisson-hurdle.
## 2 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=2000.
## 
##        mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
## lambda 4.77    0.01 0.28 4.22 4.58 4.76 4.95  5.33  1449    1
## theta  0.17    0.00 0.04 0.10 0.14 0.17 0.20  0.27  1632    1
## 
## Samples were drawn using NUTS(diag_e) at Fri Jul 12 12:30:17 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Here we compare the posterior densities of \(\lambda\) from the simple Poisson model and the Poisson hurdle model. Please finish the code:

# Extract the sampled vlaues for lambda, and store them in a variable called lambda_draws2:
lambda_draws2 <- as.matrix(fit2, pars = "lambda")

# Compare
lambdas <- cbind(lambda_fit1 = lambda_draws[, 1],
                 lambda_fit2 = lambda_draws2[, 1])

# Shade 90% interval
mcmc_areas(lambdas, prob = 0.9) 

Obtaining posterior samples from the hurdle model is more complicated, so the Stan file includes code to draw from the posterior predictive distribution for you. We extract them here, and store the predictive samples in `y_rep2’.

y_rep2 <- as.matrix(fit2, pars = "y_rep")

Plot the remaining PPCs, and comment on how this second model compares to both the observed data and to the simple Poisson model.

ppc_hist(y, y_rep2[1:8, ])

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ppc_dens_overlay(y, y_rep2[1:60, ])

ppc_stat(y, y_rep2, stat = "prop_zero")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ppc_stat_2d(y, y_rep2, stat = c("mean", "var"))

ppc_error_hist(y, y_rep2[1:4, ], binwidth = 1) + xlim(-15, 15)

## Warning: Removed 8 rows containing missing values (geom_bar).

Leave-one-out cross-validation

Here, we assess the predictive performance of both models with leave-one-out cross-validation (LOOCV). LOOCV holds out a single observation form the sample, fits the model on the remaining observations, uses the held out data point as validation, and calculates the prediction error. This process is repeated \(n\) times, such that each observation has been held out exactly once. We make use of the hand loo() function, and compare the two models:

log_lik1 <- extract_log_lik(fit, merge_chains = FALSE)
r_eff1 <- relative_eff(exp(log_lik1)) 
(loo1 <- loo(log_lik1, r_eff = r_eff1))

## 
## Computed from 2000 by 73 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -203.1 14.7
## p_loo         2.5  0.5
## looic       406.2 29.3
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

log_lik2 <- extract_log_lik(fit2, merge_chains = FALSE)
r_eff2 <- relative_eff(exp(log_lik2)) 
(loo2 <- loo(log_lik2, r_eff = r_eff2))

## 
## Computed from 2000 by 73 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -183.3 10.8
## p_loo         2.8  0.5
## looic       366.6 21.7
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

compare(loo1, loo2)

## elpd_diff        se 
##      19.8       8.5

Which model performs better in terms of prediction?

Wrapping up

Some final questions to consider:

1. Why are PPCs important?
2. Was the second model a good fit for the data? Why or why not?
3. If someone reported a single LOOCV error to you, would that be useful? Why or why not?

*Adapted from tutorial